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Abstract 


This paper proposes a non-linear state-space dynamic model for planar proton exchange membrane fuel cell (PEMFC). Our objective is the 
realization of a model that evokes a more realistic approach of dynamic behavior of the fuel cell by considering most of elements that influence 
the system evolution. For instance, pressure, temperature and humidification rate effects on the cell resistance are taken into account. The model 
is based on both thermodynamic and electrical aspects, proposing a realistic equivalent circuit which integrates most of the fuel cell components. 
Simulation results show that our proposed model is in agreement with fuel cell real operating principles. 


© 2006 Elsevier B.V. All rights reserved. 
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1. Introduction 


Since the decrease of the oil reserves, due to the massive 
exploitation of the fuel, the domain of energy becomes the 
focus of public interest. In addition, the environmental pollu- 
tion requires the development of new technologies of energy. 
Given these facts, the fuel cell becomes a prime candidate as an 
energy source and a natural alternative, which gives innovating 
answers concerning the energy effectiveness. 

In this paper, we study the proton exchange membrane fuel 
cell (PEMFC) which is considered as the most appropriate type 
of fuel cell for replacing internal combustion engines. It uses 
a solid polymer membrane (originally the Nafion 117) as elec- 
trolyte and usually works at average temperature (in the range 
of 40-90°C). Compared to other types of fuel cells, PEMFC 
is more compact and lightweight because it generates a bet- 
ter volumic-power and power-density. In addition, its operating 
temperature is less than 100°C, allowing rapid start-up. These 
facts and the ability to rapidly change power output are some of 
characteristics that make PEMFC suitable for automotive power 
applications. Other advantages result from the solid nature of the 
electrolyte because with a solid electrolyte the sealing of anode 
and cathode gases is far easier, and therefore, less expensive to 
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manufacture. It has also less problems with corrosion, which 
means a longer cell and stack life. 

There are several studies in the literature related to the 
PEMFC modeling. These studies aim at understanding the phe- 
nomena occurring in the fuel cell in order to make the appropriate 
control. A number of mathematical models have been developed 
for this purpose but they are mostly steady-state models which 
are typically used for component sizing [1,2], thermodynamical 
aspect [3-6] and electrical models [7,8]. These models repre- 
sent each component such as compressor, heat exchanger and 
fuel cell stack voltage as a static performance or efficiency map. 
Another interesting model is the one developed by [9]. Based 
on the linearized Nernst potential, the authors have presented a 
simple basic isothermal cell model which is independent of the 
type of fuel cell and is used for a general description of cells 
with gaseous fuel. 

While the above studies have been stated only for steady- 
state modeling, equivalent work have been presented for the 
dynamical modeling. Indeed, [10] uses the state-space linear 
form while [11] uses the Bond Graph method to model each 
component alone. 

In summary, a large amount of literature has been published 
about modeling of PEMFCs but most of them are centered on 
either static electrochemical or linear dynamic aspect. However, 
non-linear dynamic characteristics of fuel cell must be modeled 
in order to ensure a better dynamic control which is suitable for 
dynamic applications such as vehicles propulsion. With the best 
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Nomenclature 
A area (m?) 
C concentration (mol m~?) 


Cal, Cgeom double layer and geometrical 
capacity (F) 

D diffusion coefficient (m? s7!) 

E Nernst voltage (V) 

F Faraday constant = 96500 C mol7! 

AG Gibbs free energy (kJ mol~!) 


i output current cell (A) 
j current, density (A cm~?) 
J diffusion flux (umol s7! m7?) 
J consumption rate at triple phase boundary 
(umol s7 1) 
Js consumption rate flowing into the outer surface of 


the diffusion layer (mol s~!) 
thickness (m) 

molar mass (g mol!) 

number of exchanged electrons 
partial and total pressure (Pa) 
pore radius (cm) 

resistance (Q) 

gas constant = 8.314 J (mol K)"! 
Laplace operator 

temperature (K) 

inputs vector 

auxiliary input 

voltage (V) 

electrode voltage (V) 


a eS 2n 


(e) 


Greek letters 


a charge transfer coefficient = 1 
€ porosity 

n voltage loss (V) 

Am humidification rate 

o electrolyte conductivity (Q m)~! 
T tortuosity 

Subscripts 

a anode 

act activation 

c cathode 

conc concentration 

ct charge transfer 

load load 

m membrane 

o ohmic 

out output 

Superscripts 

b gas flow bulk 

tpb triple phase boundary 

0 at standard pressure 


of our knowledge, a notable work that takes into account this 
inherent non-linearity is the one given by [12]. In that paper, 
the authors present a cell-level non-linear state-space dynamic 
model of solid oxide fuel cell (SOFC). However, their model 
do not take into account the activation loss, the evolution of 
some parameters as well as the temperature evolution and the 
membrane humidity. They consider that electrochemical reac- 
tion potential, ohmic resistance and charge transfer resistance are 
constants while they depend respectively on pressures, humidi- 
fication rate, temperature and current. 

We propose in this paper a non-linear state-space dynamic 
model for single cell of a planar PEMFC structure that takes 
into account the influence of several parameters on the system 
behavior. Our model is based on both thermodynamic and elec- 
tric aspects. On the one hand, the diffusion principle, which is 
drawn from [12], is used to describe pressures and flow rates 
evolutions. On the other hand, the electrical principle is used 
to model the voltage and the current. The dynamic behaviors 
of voltage, current and gas consumption rates are stated with 
respect to load, partial pressures, humidification rate and tem- 
perature. Furthermore, a more realistic equivalent circuit which 
integrates most of the fuel cell components is also proposed. 

The remainder of this paper is organized as follows: Section 
2 explains the PEMFC principle and states some known results 
and remarks, which are needed in the subsequent development. 
Section 3 presents a way to build the dynamic model of a sin- 
gle PEM-type cell. Finally, simulation results and analysis are 
discussed in Section 4. 


2. Fuel cell principle 


In this section, we describe the PEMFC principle from 
reaction equation and we point to the electrochemical equations 
as well as the voltage losses and current evolution. Indeed, a 
fuel cell is an electrochemical energy conversion device which 
produces electricity, water and heat using fuel and oxygen in 
the air. It is interesting to notice that water is the only emission 
when hydrogen is the fuel. The principle of PEMFC is shown in 
Fig. 1. 

At the anode, the hydrogen molecule is split into hydrogen 
ions (protons) and electrons. The hydrogen ions permeate across 
the electrolyte to the cathode while the electrons flow through 
an external circuit and produce electric power. Oxygen, usually 
in the form of air, is supplied to the cathode and combines with 
electrons and hydrogen protons to produce water. Electrons cir- 
culation between anode and cathode through an external circuit 
generates electricity. The reaction that takes place at the triple 
phase boundary zone (TPB), is described as follows: 


H2 + 402 — H2O + electricity + heat. (1) 


At standard pressure, each element of the fuel cell (a single 
cell) makes the direct transformation of chemical energy (Gibbs 
free energy) to electrical energy according to the following equa- 
tion: 


AG +nFE® = 0. (2) 
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Fig. 1. Principle of PEMFC. 


with n = 2 in this case, because the reaction involves two elec- 
trons. 

Output voltage of the fuel cell is an important part of its 
modeling. At normal conditions of pressure and temperature, 
the electrochemical reaction potential is given by the Nernst 
equation [13]: 


E= E+ H n (2e), GB) 


PH20 


where H = RT /nF and E? = 1.273 — 2.7645 x 10~4T. 
When a fuel cell operates, the actual voltage is lower than 
the one computed with Nernst’s equation, due to activation, 
diffusion and ohmic losses. The reduction of the voltage is pro- 
portional to the current circulation through the cell. Fig. 2 shows 
the voltage/current characteristic taking into account losses [14]. 
Fig. 2 shows three evolution zones for the output voltage: 


e The first zone represents the activation loss caused by the 
charge transfer between electrodes and electrolyte. It appears 
only for a little current and is governed by the following 
equation: 


Vout (V) 


J (A/em2) 


Fig. 2. Polarization curve. 


RT j 
Nact = hn >], (4) 
anF Jo 


where jo denotes exchange current density. 
e The second zone represents the ohmic loss due to the resis- 


tance of different components of the fuel cell. Its expression 
is as follows: 


No = Rott, (5) 


where i; is the total current. 

e The third zone represents the concentration loss due to the 
mass transfer of reagents. It appears only for a huge current 
and is governed by the following equation: 


RT j 
Neone = —— Inj 1- >], (6) 
anF j 


where jı denotes limit current density. 
Hence, the effective output voltage of the fuel cell is given 
by the following equation: 


Vout = E — nact — No — Neone- (7) 


Now, concerning the current output, the thermodynamic 
aspect of the fuel cell allows us to establish the relation between 
current and gas flow rates as follows [15]: 


i = 2F JG, = 2FJipo = 4FJ6,. (8) 


It is interesting to note that the current is not only influenced 
by load value but it depends on another variables. This fact can 
be easily seen in Eq. (8) which shows that current depends on gas 
flow rate and thus on gas pressures values. Besides, the current is 
limited by the maximum flow rate value which in turn depends 
on concentration gradients between TPB and gas flow bulks. 
When the current output increases, the hydrogen and oxygen 
concentrations at TPB decrease to create larger concentration 
gradients. However, because the cathode reaction is slower than 
the anode one, current output is also limited by the maximum 
ion production rate. In addition, the current is influenced by the 
membrane resistance which depends on the temperature and the 
humidification rate. 


3. Modeling approach 


In general, two types of fuel cell models have been considered 
in the literature: steady-state and dynamic models. The former 
focuses on simulating the fuel cell polarization curve while the 
latter generally pays more attention to thermodynamic aspects. 
This paper deals with modeling the dynamic performances of a 
fuel cell focusing on both electrical and diffusion approaches. 
The dynamics of the PEMFC will be represented by non-linear 
differential equations transformed into an equivalent mathemat- 
ical formulation: the state-space representation. The method of 
state-space domain exhibits great advantages especially in the 
study of control, parameter estimation, data reconciliation and 
optimization. Besides, it presents a powerful tool for simulation. 

As mentioned in Section 1, we are interested on the planar 
fuel cell geometry. This structure is built by stacking together 
multiple single cell. Hence, the fuel cell stack is electrically a 


A. Haddad et al. / Journal of Power Sources 163 (2006) 420-432 423 


Table 1 
Input/output variables 
Inputs Outputs 
u1 = Road Load resistance yı = Vout Voltage output 
u = pe, Partial pressure of hydrogen y2 =i Current output 
in anode gas bulk 
u3 = p>, Partial pressure of oxygenin y3 = Jis Hydrogen 
cathode gas bulk consumption rate 
u4 = Po Partial pressure of water y= Jb, Oxygen 
vapor in anode gas bulk consumption rate 
ug = À Humidification rate ys = Jħįo Water 
consumption rate 
us =T Fuel cell temperature 


serial connection of these cells. It is basically this fact which 
motivates us to focus on a single cell modeling since in a serial 
connection the single cell voltage is added up. In addition, we 
consider that the fuel cell lies in a mesoscopic scale. For this 
reason, we may consider a single cell as a lumped volume type 
OD of the stack. This enables us to assume an uniform variation 
of the temperature inside the single cell. However, it is important 
to note that the study of temperature variation with respect to 
the space scale needs a special thermal study which is not the 
subject of the paper. 

According to these hypothesis, we develop a model based 
on diffusion and electrical principles. The diffusion principle 
enables us to identify pressures and flow rates while electrical 
principle allows us to model the voltage and the current. The 
PEMFC modeling process will be divided into three parts. In 
the first part, we establish relations between pressures and flow 
rates according to concentration and diffusion laws. Then, we set 
electrical relationships between voltage and current in the sec- 
ond part. The third part is devoted to the state-space formulation. 
To do this, we define in Table 1 the input/output variables. 

The first input is the external load that influences the current 
output i and so affects the reaction. In order to investigate the 
basic dynamic behavior of the current output, load impedance is 
assumed to be a pure resistance Rjoaq in the model. Other inputs 
are partial pressures of reactants in gas bulks. They affect the 
diffusion processes as well as the Gibbs free energy and thus the 
voltage and current outputs. The fifth input is the humidification 
rate A related to the water management. Its role is to maintain the 
membrane humidity and to balance water usage/consumption in 
the system. Indeed, on one hand, much water leads to the flood- 
ing of the membrane, which in turn blocks the mass diffusion 
and proton transportation. On the other hand, if the humidifica- 
tion level of the membrane is too low, the proton transportation 
through the membrane will become difficult and this will lead 
to higher internal resistance. Another input that influences the 
membrane resistance is the temperature T which together with 
humidification rate influences the ohmic losses. This eventually 
affects voltage and current outputs. Finally, the output variables 
are voltage, current and gases consumption rates. 


3.1. Diffusional approach 


In this part, we are going to use the diffusion principle in 
order to establish the expressions of gas partial pressures p'P? 


at TPB and gas consumption rate J at gas bulk. As the gas 
consumption rate J’ at TPB can be expressed with respect to 
the current (see Eq. (8)) and p? is taken as an input variable, we 
will establish the expressions of p'P® and J’ with respect to both 
p? and J". Indeed, in our planar geometrical model for which 
the co-flow circulation of gases is considered, we determine the 
diffusion flux due to the concentration difference between gas 
bulk and TPB zones according to Fick’s first law: 


J = —DVC, (9) 


where VC is the concentration gradient. The low diffusion layer 
thickness and the porous characteristic of the electrodes allow us 
to assume that gases diffuse in one dimension (in the x direction; 
see Fig. 1). Hence, the expression of the diffusion flux Eq. (9) 
becomes: 


J =—-D—. (10) 


The diffusion coefficient D in porous materials can be calcu- 
lated according to Giavazzi—Pagano simplified equation: 


€e IT 
D = 9700r -4| — 
TV M 


In order to compute J we use Fick’s second law which pro- 
vides us the variation of the concentration with respect to time: 


ac 3C 
— = D—. 

ot ax2 
For the solution of this equation, as in [12], we use the Laplace 
transforms which appears to be more suitable for the state-space 


representation. Indeed, the Laplace transforms of Eq. (11) is 
given by 


(1) 


Chs) s Cis) <0 (12) 
dx? D peor Ss 
with the following boundary conditions: 
dc 
ronan 
dx x=0 
dC 
fee ap (13) 
dx x=L 


C(s) = C()lxaz 


Hence, the solution of Eq. (12) with respect to these boundary 
conditions is as follows: 


J™(s) sinh(./s/Dx) 
AD\/s/D 
J"(s) sinh(./s/DL) + C°(s)AD./s/D 
AD,/s/D cosh(./s/DL) 


x cosh (5) A (14) 


Now, assume that gases are ideal. Then we have: 


NRT 


C(s)\(x) = 


PV = NRT => P = = CRT, (15) 
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where N is the number of moles of gas present on the volume 
VY. At TPB we have x = 0 (see Fig. 1). This infers that p'P>(s) = 
RTC(s)|x=0. With this in mind, from Eq. (14) we derive: 


Cts) = sinh(ys/DL) _, Ow 
om = ADJs/D enue | cosh(/s/ DL) 


This and the fact that p?(s) = RTC(s)|,—o imply: 
RT sinh(./s/DL) 
AD,/s/D | cosh(./s/DL) 
p(s) 
cosh(./s/DL) 


Now, at the gas bulk zone we have x = L. It follows after a 
simple calculation that 


p”™(s) = 


| J™(s) 


(16) 


Bice: J"(s) AD4/s/D Vs ) b 
an cosh(./s/DL) a RT an ( D Bits). 


(17) 


As in [12], since L is around zero, the use of the Taylor’s 
serie expansions of cosh, sinh and tanh functions about 0 and 
neglecting higher order terms we arrive at: 


p(s) = G yp J*(s) + G pp p(s), 


‘ (18) 
J8(s) = Gy J™(s) + Gps p(s), 
where 
a —(L/D) —(L3/6D?)s RT, 
T+ (L2/2D)s + (L4/24D2)s2 A’ 
1 
Om= ig (L2/2D)s + (L4/24D2)s2 
1 
Gyy= ; 
IT T+ (L2/2D)s + (L4/24D2)s2 
Ls A 


Gp; = 
PI T+ (L2/2D)s + (L4/24D2)s2 RT 
3.2. Electrical approach 


It is by now well known that representing the fuel cell sys- 
tem as an electrical circuit is quite useful and a great effort has 
been made to assimilate the fuel cell to an inherent impedance 
by modeling the electrical effects of its components. Accord- 
ing to the PEMFC principle (Section 2), fuel cell’s voltage 
output is affected by gas partial pressures and is reduced by 
activation, concentration and ohmic losses. Because reactions 
take place at TPB, the more appropriate expression of Nernst 
equation is: 


(19) 


The electrochemical phenomena produced inside the fuel cell 
components are modeled by electric impedances. The inher- 
ent impedance is the equivalent circuit of these components. It 
represents the losses and its dynamic characteristics affect the 
dynamic behaviors of voltage and current. More precisely, it is 
composed of: 


Fig. 3. Equivalent circuit of inherent impedance. 


(1) Three types of resistances representing the ohmic resis- 
tances (Ro, Ra, Ro) the activation resistances (Racta, Ractc) 
and the charge transfer resistances (Reta, Retc). 

(2) Two charge capacities (Cala, Caic) due to the double layer 
phenomenon and one geometrical capacity Cgeom between 
the electrodes. Indeed, the structure of charge accumulation 
and the charge separation always occur at the interface when 
an electrode is immersed into an electrolyte solution. The 
excess charge on the electrode surface is compensated by 
an accumulation of excess ions of the opposite charge in 
the solution. This structure behaves essentially as capaci- 
tor and leads to (Cala, Caic). For Cgeom, it follows from the 
fact that two electrodes (anode, cathode) separated with an 
isolator (electrolyte) produce a capacitance effect between 
them. 


In fact, the inherent impedance Zinherent is modeled as an 
equivalent RC circuit. A typical equivalent circuit which takes 
into account all inherent impedances is shown in Fig. 3. 

Introducing the electrochemical reaction potential and the 
load resistance to the inherent impedance we obtain the equiv- 
alent circuit of the PEMFC shown in Fig. 4, where i, is the 
current of the geometrical capacity. The model shown in Fig. 4 


Fig. 4. Equivalent circuit of one cell. 
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is an integral circuit that appears to be realistic by introducing 
activation, ohmic and concentration losses as well as the double 
layer and geometrical capacitance effects. Resistances are calcu- 
lated with respect to the inputs while capacities can be identified 
from impedance spectral plan plot [11]. 

Let Ve be either V, or Ve. According to Ohm’s law, the expres- 
sion of Ve from the equivalent circuit is as follows: 


dV. 
= —(Ract + Ret) | it + Ca- l (20) 
This implies that: 
dVe —1 1 
= Ve 
dt (Ract + Ret)Cai Cai 


Therefore, the voltage and current outputs are calculated as 
follows: 


it. (21) 


Vet =E+ 2Ve 
Vout = Vet = itRo =E+ 2Ve = itRo (22) 
aa 4 i Vout LC dVout 
t= C= ima geom `~ | dt 
This leads to: 
Ve dV, 
Vout = E + 2Ve — Ro a + Cgeom—— out (23) 
load dt 
The solution of Eq. (23) with respect to the initial condition 
Vout(0) = 0 gives: 
Rioad(E + 2Ve)[1 — exp(—((Rioaa + Ro)/ 
Rioad RoC t 
Vout(t) = load No geom) )] (24) 


Rioad + Ro 
3.3. State-space representation 


In this part, we present the model based on diffusional and 
electrical approaches in the state-space form. In order to model 
the dynamic behavior of gas consumption rates, we convert the 
transfer function form the dynamic relations shown in Eq. (18) 
to a differential equation form as follows: 


e Hydrogen consumption rate: 


Ji, = PRL 


78 r 
4 —a Ji, — a2 Sq, + a1 Sq, + 03 


RT 


and its partial pressure in the vicinity of anode at TPB: 


ut tpb tpb RT 
Pu, = 1 Ph, a 2PH, a a Ju, 
4 RT. 
g T 
— — J 
La A H, + 1PH, 


where a = 24D, /L4, a2 = 12Dy,/L3, 03 = 24D, /L3, 
a4 = 24Dp, / L3. 
e Oxygen consumption rate: 


"7S N TS T A -b 
Jo, = —b1 Jò, — b2Jò, + B1 Jo, + P3 p r ËO 
8 


and its partial pressure in the vicinity of cathode at TPB: 


tpb tpb tpb ReT 
BO, = —B1Pby — BrP; — Ba To, 
4 Rely +B 
TR A 1PÒ, 


where 6; = 24D, /Lé, Bo = 12Do,/L?2, B3 = 24D, ji. 
ba = 24Do,/L?. 

e Water vapor production rate: 

—yı Jibo = 


7S TS T A, b 
Jio = v2Jņ o + Y1 Jm o + v3 R TPO 
g 


and its partial pressure in the vicinity of anode at TPB: 


tpb Rg T r 


Po = = Pino ~ ¥2PH50 — y4— y o 


4 RT., i 
= L, A #20 + yı Pho 
where yı = 24Di,,9/L3, ¥2 = 12Dyn0/L}, y3 = 24Di,0/ 


L, y4 = 24Dy,0/L?3. 


The above equations show that gases consumption rates 
depend on input derivatives (ù2, ù3 and w4). This fact is gen- 
erally negative because the derivative mode can amplify the 
outputs noise. Therefore, we use a low-pass filter to eliminate 
the noise increasing in derivative mode. Consequently, the first 
order derivative of the input variables can be approximated by 
the following equation [16]: 


1 
sU(s) ~ K (1 ete z) U(s), 


where K is an approximation factor greater than 10. Let v be an 
auxiliary input defined as v(s) = (1/(1 + s/K))U(s). Then we 
have: 
ù = Ku — v 
ù = K?u — Kv 

Hence, after introducing the intermediate variables vp, , vo,, 


UH»O and vR,,,, to the model, we define the state-space vector as 
follows: 


= tpb tpb tpb tpb ptpb pie 
x = [Vex PH» PO Pho: Py >» PO, + Pibo Jio Stig UH» JO,» 


7S s s T 
Jò,» voz; Jino: Jio» VH20; VRioad] 


This together with all the above equations lead to the follow- 
ing state-space model: 


x) it 


a (Ract + Ret)Car Ca 

ko = x5; 

x3 = X63 

x4 = X7; 

A Regie . 2Rgu6 dit 
X5 = —Q1X2 + &1U2 — 2X5 — a4 ZFA lt FAL, ar? 
že = —Bix3 + Biu3 — Boxe — p4 sei ae z 
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. 7 Regio. 2Reu dit 
x7 = —y1x. u x i : 
a a Se ay a FAL, di 

xg = x9; 

ial mca e « ) 

Xo = —Q1 Xg — ADX l u X10); 

9 1X8 2X9 F alt Relg 2 10 

x10 — Ku = Kx103 

X11 = x12; 

; Bi. A 

X12 = — 1x11 — 2x12 + lt + (Kuz — x13); 
4F Rete 

B= K?u3 — Kx13; 

X14 = x15; 

À yi. nA 

x15 = —yY1Xx14 — Yox15 + tt t (Ku4 — x16); 
2F Rete 


16 = K7u4 — Kx16; 
*17 = Ku, — Kxı7. 


The output equations are represented as follows: 


ui(E + 2x1)[1 — exp(—((u1 + Ro)/u1 RoCgeom)t)] | 
uj + Ro i 


y4 = x11; 


y= 


yı 
y=; 
uy 


y3 = X8; y5 = X14. 


Here, the variables i;, di,/dt, Ract, Ret, E and Ro are as fol- 
lows: 


pE (E + 2x1)[1 — exp(—((u1 + Ro)/u1 RoCgeom)t)] 


ui + Ro 
$ (E + 2x) )lexp(—(u1 + Ro)/u1 RoCgeom)t)] 
Ro 
dit = (E + 2x; )[exp(—(u1 + Ro)/u1 RoCgeom)t)] 
dt u1 RoCgeom 


(E + 2x1)(u1 + Ro)lexp(—((u1 + Ro)/u1 RoCgeom)t)] 
u1 R2Cgeom 


o = (5139 x 10-°us — 326 x 1075) e!?68(0/303)=(1/u6)), 


Hence, the state-space model is developed in the following 
form: 


x= f(x, u,t) 
y = g(x, u, t) 


This form is a non-linear dynamic model which is essential 
for the optimal operation and the control of many processes. It 
presents a powerful tool for simulating the real behavior of the 
system because the state-space form allows a detailed analysis 
of the cell temporal and spatial behavior. The temporal behavior 
provides the evolution of all state variables with respect to the 
time axis while the spatial behavior provides a so-called phase 
space which consists on the representation of some or all state 
variables in a multidimensional space. This phase space may 
give qualitative information about the dynamics of the system. 
In addition, the state-space form is suitable for the optimiza- 
tion of physical parameters and is well suited for synthesis of 
a controller to stabilize the cell voltage or reject the possible 
environmental disturbances. 


4. Simulation results 


In this section, we will present the simulation results that 
demonstrate the effect of the load change and its disturbance, the 
gas pressures change as well as the effects of humidification and 
temperature. The simulation is done according to the developed 
model and it investigates steady-state and transient behaviors of 
a single cell at different inputs and disturbances. It is done under 
Matlab-Simulink using the well known S-function toolbox as 
shown in Fig. 5. 

Since data related to dynamic models are limited in the lit- 
erature, we will simulate the model according to the parameters 
of [11,12]. Most of these data are taken from [11] because the 
author studies the PEMFC. However, some other parameters 
are taken from [12] like the cell dimension (A, La, Le) and the 
porous characteristics of the electrodes (e, t). The simulation 
results represent the dynamic behavior for one cell of PEMFC. 
The parameter values are shown in Table 2. 


4.1. Effect of load 


The load is usually imposed by the external circuit. How- 
ever, we consider it as an input in order to view the effect 
of its change on the outputs. We will present its effect by 
simulating the responses to pulse, triangular and sinusoidal 
resistance change as well as the effect of the disturbance 
(Figs. 6-11). 

One can immediately observe that the variation of the cur- 
rent is opposite to the variation of the load resistance while the 
output voltage change is proportional to it. This is in agreement 
with Eq. (22). In addition, we notice that the current variation 
is bigger comparing to the voltage one, this is due to the fact 
that the influence of the load on the current is larger than its 
influence on the voltage. Indeed, from Eqs. (22) and (24) we 
have: 


i- Vout 
Rioad i 
Rioaa( E + 2Vo)[1 — exp(—((Rioaa + Ro)/ (25) 
Vout(t) = Rioad RoC geom)t)] 
= E Rioad + Ro 
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Fig. 5. Simulink model. 

Table 2 
Parameters value 
Symbol Description Symbol Description 
A = 1cm? Fuel cell effective area Pu, = 2.5 x 10° Pa Hydrogen partial pressure in anode gas bulk 
Ca = 186 mF Double layer capacitance Po, = 2.5 x 10° Pa Oxygen partial pressure in cathode gas bulk 
Cgeom = 186 mF Geometrical capacitance Phyo = 1 x 10°Pa Water partial pressure in cathode gas bulk 
e=0.4 Porosity R, = 45mQ Anode resistance 
io = 0.2A Activation current R: = 45mQ Cathode resistance 
i =2A Limit current r= 20pm Pore radius 
La = 1mm Thickness of anode diffusion layer tT=4 Tortuosity 
L. = 1mm Thickness of cathode diffusion layer 
Lm = 0.2mm Thickness of membrane 


Because Ro < Rioaq We have Ro + Rioad © Rioaa. Given this 
fact, the expression of Vout becomes: 


Rioad(E + 2Ve)[1 — exp(—(Rioad/ Rioad RoC geom)t)] 
Rioad 


t 
Ro Ceeom ) | l (26) 


Vout(t) ad 


= (E + 2Ve) i exp ( 


— Rload (Ohm) 
== Vout (V) 


t(s) 


Fig. 6. Responses to pulse load changes. 


Hence, one can immediately observe from Eqs. (25) and (26) 
that the load influence on the current is more important than its 
influence on the voltage. 

Furthermore, simulations in Fig. 6 show that the time 
response between the input change and the output responses 
is about 0.25 s (see Fig. 9). Comparing this result with one pre- 
sented in [12], where the time response is 0.08 s, we see that 
the difference is due to the fact that in our study we take into 


— Rload (Ohm) 
== Vout (V) 


t(s) 


Fig. 7. Responses to triangular load changes. 
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— Rload (Ohm) 
-— Vout (V) 


t(s) 


Fig. 8. Responses to sinusoidal load changes. 


== Vout (V) 


1.14 ye 


t (s) 


Fig. 9. Time response of the outputs. 


account the effect of the geometrical capacity while the work of 
[12] do not consider it. 

As in [12,11], we do not have any overflow in the output 
responses. This is related to the input values. In general, the 
overflow is proportional to the velocity of the load variation. 


4F T T T 4 
3.5/5 4 
3r 4 
2.57 | 
2r 
1.5 4 
1F -d 
0.5F 4 
1 1 1 A 
0 2 4 6 8 10 


t(s) 


Fig. 10. Load disturbance. 


t (s) 


Fig. 11. Responses to load disturbance. 


— Pressure (Pa) 
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Fig. 12. Gase pressure change. 


The results of Fig. 11 show the effect of the load disturbance 
(Fig. 10) on current and voltage. 

One can observe that variables are immediately influenced by 
load disturbance, but the current is more disturbed than the volt- 


t(s) 


Fig. 13. Responses to hydrogen pressure pulse. 
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Fig. 14. Responses to oxygen pressure pulse. 
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Fig. 17. Membrane humidification rate change. 
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Fig. 15. Responses to water pressure pulse. 


age because the load influence on the former is more important 
than its influence on the latter (see explanation in this section). 
For a load disturbance of 100%, the current change is about 50% 
while the voltage one is about 4%. 


t(s) 


Fig. 16. Effect of membrane humidity. 


t(s) 


Fig. 18. Effect of temperature. 


4.2. Effect of partial pressures 
Usually, the control of the fuel cell system is not based on 


reactant pressure action, because this parameter does not have a 
big influence on the outputs. This fact will be shown in this part 
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Fig. 19. Temperature change. 
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Fig. 20. Response to Rioaa steps from 4 to 2 Q. 
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Fig. 21. Response to py, steps from 0.97 to 1.2 atm. 
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Fig. 22. Effect of temperature. 


where we simulate the effect of hydrogen, oxygen and water 
pressures by simulating the output responses to pulse pressure 
change (Fig. 12). 

According to Nernst equation, we see in Figs. 13—15 that 
the voltage response is proportional to hydrogen and oxygen 
pressure evolutions but it is conversely proportional to the water 
pressure. Comparing Figs. 13—15, we observe that the variation 
of the voltage output with respect to Peo and Ph, is more 


important than its variation with respect to Pò, . This is coherent 
with the fuel cell principle because cathode reaction is slower 
than anode one. 


4.3. Effect of humidification 


As described before, the effect of the membrane humidity is 
very important because it influences the loss behavior. Indeed, 
the membrane humidity affects its resistance value which is the 
main element that determines the ohmic loss. The simulation 
results in Fig. 16 show the output responses to the humidifi- 
cation rate change (Fig. 17) between 50 and 100%. One can 
show that the voltage Vout and the current i behave propor- 
tionally to the humidification rate Am. This is due to the fact 
that the membrane resistance Rm is conversely proportional 
to Am. 


4.4. Effect of temperature 


The fuel cell temperature is an important factor that influences 
the outputs because it affects the Nernst voltage as well as the 
activation, ohmic and concentration resistances. The simulation 
results in Fig. 18 show the outputs responses to the temperature 
change (Fig. 19) between 40 and 90°C. 

According to Eqs. (3) and (22), we can see in Fig. 18 that both 
voltage and current are proportional to the temperature. This is 
due to the fact that Nernst voltage and membrane resistance are 
proportional to the temperature. 

In conclusion, we can emphasize the fact that the dynamic 
behavior of the fuel cell outputs is mainly influenced by the fuel 
cell temperature, the humidification of the membrane and the 
external load. 


5. Comparison with other work 


As we mentioned in the introduction, most of the work in the 
literature are centered on either static electrochemical aspect 
or linear dynamic aspect. With the best of our knowledge, the 
notable work which takes into account the non-linearity aspect is 
the one presented in [12]. However, the paper focuses on the solid 
oxide fuel cell and do not take into account the activation loss, 
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the geometrical capacity as well as the evolution of some impor- 
tant variables like the fuel cell temperature and the membrane 
humidity. In addition, we do not have the same functional condi- 
tions, for instance, the functional temperature for the PEMFC is 
80°C while it is 800 °C for the SOFC. For all these reasons, our 
output values are different comparing to those presented in [12] 
but nevertheless they must have the same behavior due to the 
fuel cell principle. The comparison diagrams show the behavior 
similarity of Vout and i with respect to Rjoad, PH, and T change 
in both cases. The only important point that should be noted is 
that in our results (Figs. 20-22, left side) the overflow appears 
at the pressure response contrarily to the results of [12] (Figs. 
20-22, right side) which appears at Rjoad response. Besides, we 
do not consider the same time scale because we have a bigger 
time response. 


6. Conclusion 


In this paper, a non-linear dynamical state-space model of 
planar proton exchange membrane fuel cell has been built. In 
our model, we have considered the majority of the phenomena 
that influence the dynamic behavior of the fuel cell. In partic- 
ular, a realistic equivalent circuit has been proposed. In that 
circuit, we have considered all the fuel cell elements that cause 
the energy losses. Furthermore, the pressure effect on the electro- 
chemical reaction potential as well as the one of the temperature 
and the humidification rate on resistances have been also taken 
into account. Simulation results show that the main elements 
that influence the outputs are temperature, humidification rate 
and load. Hence, our model is in agreement with the fuel cell 
principle. 

In conclusion, with the aid of dynamical model, there are 
several issues that deserve further investigation. One is the con- 
trol of the energy losses in the fuel cell. Indeed, among the real 
obstacles (such as production costs and life limitations) for fuel 
cells to become a success is their energy losses. This problem 
is even more accentuated for automotive applications because 
the energy loss lies at the core of performance problems. The 
performance of a fuel cell depends not only on the individual 
components but also on how the components are controlled. 


Hence, minimizing energy losses can be applied to the model 
using methods from optimal control. This is a topic of our future 
research. 
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